We introduce time steps
\( t_0=0 < t_1 < t_2 < \cdots t_{k-1} < t < t_{k+1}< \cdots < t_N \)
and states of the modeled phenomenon in particular time steps, represented by a scalar field \( u(t) \)
\( u_0=u(t_0), u_1=u(t_1), u_2=u(t_2), \cdots, u_{k-1}=u(t_{k-1}), u_k=u(t_k), u_{k+1}=u(t_{k+1}), \cdots, u_N=u(t_N) \).
We approximate the derivative in time using the finite difference method \( \frac{\partial u}{\partial t} \approx \frac{u_{t+1}-u_t}{dt } \).
So we get the equation
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}(u) = f \).
The implicit method, called the Crank-Nicolson method (invented by Mrs. Phyllis Nicolson, a female British mathematician and John Crank, an American mathematician and physicist in 1947) works as follows [1].
In this method, we apply our differential operator to the mean value of the previous and next states:
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}(\frac{1}{2}(u_t+u_{t+1})) = {f_{t+1}} \).
If the differential operator is linear, then we can break it into two terms \( u_{t+1} -dt \frac{1}{2}\mathcal{L}(u_{t+1}) = u_t+dt{f_{t+1}}+ dt\frac{1}{2}\mathcal{L}(u_t) \)
leaving our "unknown" on the left-hand side \( u_{t+1} \).
Then we use the isogeometric finite element method.
To develop a solver that uses the method implicitly in the isogeometric finite element method, we need to transform a strong formulation into a weak formulation.
So we multiply our equation by the test functions \( (u_{t+1},w) -dt \frac{1}{2}(\mathcal{L}(u_{t+1}),w) = (u_t,w)+dt({f_{t+1}},w)+ dt\frac{1}{2}(\mathcal{L}(u_t),w) \)
We will use a linear combination of the B-spline function to approximate the state of our system at a given moment in time. For this purpose, we select the basis of two-dimensional B-spline functions, defining the node vectors in the direction of the coordinate system axis.
For example, we can choose a two-dimensional basis of a second order B spline
\( \{B^x_{i,2}(x)B^y_{j,2}(y)\}_{i=1,...,N_x;j=1,...,N_y} u_{i,j } \)
These will be used to approximate the simulated scalar field to the current time instant
\( u(x,y;t+1)\approx \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } B^x_i(x)B^y_j(y) \)
Similarly, for testing, we will use the B-spline base functions:
\( w(x,y) \in \{ B^x_k(x)B^y_l(y) \}_{k=1,...,N_x;l=1,...,N_y} w^{k,l } \)
Assuming that the differential operator describing "physics" is linear, that is \( {\cal L}(\alpha u+\beta w)=\alpha {\cal L}(u)+\beta {\cal L}(w) \) our equation now looks like this:
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } (B^x_i(x)B^y_j(y)-dt \frac{1}{2} \mathcal{L}(B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y))=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \)
For example, for the problem of heat transport we have
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ (B^x_i(x)B^y_j(y)-dt \frac{1}{2} \left(\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2 }\right),B^x_k(x)B^y_l(y))=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \)
Thanks to the weak formulation, we can integrate by parts
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ \left( (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) - dt \frac{1}{2}(\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x},\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}) -dt* (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y }) \right)=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \)
Due to the structure of the tensor product of the B-spline function and due to the fact that the derivative in the y direction of the B-spline in the x-axis direction gives 0 (because these functions are constant in the y-axis direction) and similarly for the derivative in the y-direction of the B-spline in the x-axis direction, we have
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ \left( (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) - dt \frac{1}{2}(\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y)) -dt* (B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y }) \right)=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \)
Note that our system of equations to solve is \( \left(M_x \otimes M_y-dt \frac{1}{2} S_x \otimes M_y -dt* M_x \otimes S_y\right) u^{t+1} = F(t) \)
where
\( \{ M_x \otimes M_y \}_{i,j,k,l}= \int B^x_i(x)B^x_k(x)B^y_j(y)B^y_l(y) = \int B^x_i(x)B^y_j(y)B^x_k(x)B^y_l(y) = {\bf M}_{i,j,k,l } \) is a mass matrix which is a Kronecker product of two one-dimensional mass matrices,
\( \{ S_x \otimes M_y\}_{i,j,k,l} = \int \frac{\partial B^x_i(x)}{\partial x}\frac{\partial B^x_k(x)}{\partial x}B^y_j(y)B^y_l(y) = \int \frac{\partial B^x_i(x)}{\partial x}B^y_j(y)\frac{\partial B^x_k(x)}{\partial x }B^y_l(y) \) is the Kronecker product of the one-dimensional stiffness matrix and the one-dimensional mass matrix, and
\( \{ M_x \otimes S_y \}_{i,j,k,l} = \int B^x_i(x) B^x_k(x) \frac{\partial B^y_j(y)}{\partial y} \frac{\partial B^y_l(y)}{\partial y } = \int B^x_i(x)\frac{\partial B^y_j(y)}{\partial y}B^x_k(x)\frac{\partial B^y_l(y)}{\partial y } \) is the Kronecker product of the one-dimensional mass matrix and the one-dimensional stiffness matrix.
Each of these matrices can be factored in linear time using the alternating-directional solver algorithm. However, it is no longer possible to factorize their sum in linear time. This is possible only when we introduce a time step scheme that allows the matrix to be split into sub-matrices in time sub-steps so that the factorization cost remains linear.
The Peaceman-Reachford scheme allows you to split the time step into two sub-steps \( \left(M_x \otimes M_y-dt \frac{1}{2} S_x \otimes M_y\right) u^{t+1/2} = F(t+1/2)+\left(dt \frac{1}{2} M_x \otimes S_y\right) u^{t} \)
\( \left(M_x \otimes M_y-dt \frac{1}{2} M_x \otimes S_y\right) u^{t+1} = F(t+1/2)+\left(dt \frac{1}{2} S_x \otimes M_y\right) u^{t+1/2 } \)
where we can merge the left-hand side matrices into a single matrix with a Kronecker product structure that can be factored in linear time using the variable direction solver algorithm
\( \left(M_x-dt \frac{1}{2} S_x \right)\otimes M_y u^{t+1/2} = F(t+1/2)+\left(dt \frac{1}{2} M_x \otimes S_y\right) u^{t} \)
\( M_x \otimes \left( M_y-dt \frac{1}{2} S_y\right) u^{t+1} = F(t+1/2)+\left(dt \frac{1}{2} S_x \otimes M_y\right) u^{t+1/2 } \)
Finally, let us take a look at what the right-hand side operator looks like.
The right-hand side function here represents the changes caused by the force acting on the system during the time step, plus the changes caused by the physics of the phenomenon in the previous time step \( (u_{t+1},w) -dt \frac{1}{2}(\mathcal{L}(u_{t+1}),w) = (u_t,B^x_k(x)B^y_l(y))+dt({f_{t+1}},B^x_k(x)B^y_l(y))+ dt\frac{1}{2}(\mathcal{L}(u_t),B^x_k(x)B^y_l(y)) \).
In particular, our right-hand side is not a bitmap, but a projection of the sum of three elements:

1. The state of our system at the previous moment in time \( (u_t,w)= \sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) \) (also multiplied by the test function and area integrated). Of course, we also use a linear combination of B-spline basis functions to represent the state of our system in the previous time step \( u(x,y;t)=u_t=\sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } B^x_i(x)B^y_j(y) \).
2. Changes induced by the "physics" of the simulated phenomenon during the time step. These changes are calculated by applying the differential operator representing the modeled phenomenon to the state of the system at the previous moment. The state of the system is of course represented by a linear combination of the B-spline function. So the differential operator is applied to the B-spline function \( (dt \frac{1}{2} \mathcal{L}(u_t),w)=\sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } ({\cal L} (B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y)) \). For example, for the problem of heat transport we have \( \sum_{i=1,...,N_x;j=1,...,N_y} dt \frac{1}{2} u^{t}_{i,j } (\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2},B^x_k(x)B^y_l(y)) \). Thanks to the weak formulation, we can integrate by parts \( \sum_{i=1,...,N_x;j=1,...,N_y} dt \frac{1}{2} u^{t}_{i,j } \left( (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x}B^y_j(y),\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},B^x_k(x)\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y})\right) \) which due to the structure of the tensor product of the B-spline function, and due to the fact that the derivative in the y direction of the B-spline function of the variable x gives 0, and the derivative in the x direction of the B-spline function of the variable y gives, we get \( \sum_{i=1,...,N_x;j=1,...,N_y} dt \frac{1}{2} u^{t}_{i,j } \left( (\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y})\right) \). From the "matrix" point of view, this term gives us \( dt \frac{1}{2} (S_x \otimes M_x +M_x\otimes S_y) u_t \).
3. Changes caused by a force acting on the system during the time step \( (f_{t+1},w)= (f_{t+1}(x,y),B^x_k(x)B^y_l(y)) \).